module clm_initializeMod

  !-----------------------------------------------------------------------
  ! Performs land model initialization
  !
  use shr_kind_mod    , only : r8 => shr_kind_r8
  use shr_sys_mod     , only : shr_sys_flush
  use shr_log_mod     , only : errMsg => shr_log_errMsg
  use spmdMod         , only : masterproc
  use decompMod       , only : bounds_type, get_proc_bounds, get_proc_clumps, get_clump_bounds
  use abortutils      , only : endrun
  use clm_varctl      , only : nsrest, nsrStartup, nsrContinue, nsrBranch
  use clm_varctl      , only : is_cold_start, is_interpolated_start
  use clm_varctl      , only : iulog
  use clm_varctl      , only : use_lch4, use_cn, use_cndv, use_c13, use_c14, use_fates
  use clm_varctl      , only : use_soil_moisture_streams
  use clm_instur      , only : wt_lunit, urban_valid, wt_nat_patch, wt_cft, fert_cft, irrig_method, wt_glc_mec, topo_glc_mec, haslake
  use perf_mod        , only : t_startf, t_stopf
  use readParamsMod   , only : readParameters
  use ncdio_pio       , only : file_desc_t
  use GridcellType    , only : grc           ! instance
  use LandunitType    , only : lun           ! instance
  use ColumnType      , only : col           ! instance
  use PatchType       , only : patch         ! instance
  use reweightMod     , only : reweight_wrapup
  use filterMod       , only : allocFilters, filter, filter_inactive_and_active
  use FatesInterfaceMod, only : set_fates_global_elements
  use dynSubgridControlMod, only: dynSubgridControl_init, get_reset_dynbal_baselines
  use SelfTestDriver, only : self_test_driver

  use clm_instMod
  use SoilMoistureStreamMod, only : PrescribedSoilMoistureInit
  ! 
  implicit none
  private  ! By default everything is private

  !
  public :: initialize1  ! Phase one initialization
  public :: initialize2  ! Phase two initialization
  !-----------------------------------------------------------------------

contains

  !-----------------------------------------------------------------------
  subroutine initialize1(dtime, gindex_ocn)
    !
    ! !DESCRIPTION:
    ! CLM initialization first phase
    !
    ! !USES:
    use clm_varpar       , only: clm_varpar_init, natpft_lb, natpft_ub, cft_lb, cft_ub, maxpatch_glcmec, nlevsoi
    use clm_varcon       , only: clm_varcon_init
    use landunit_varcon  , only: landunit_varcon_init, max_lunit
    use clm_varctl       , only: fsurdat, fatmlndfrc, noland, version  
    use pftconMod        , only: pftcon       
    use decompInitMod    , only: decompInit_lnd, decompInit_clumps, decompInit_glcp, decompInit_lnd3D
    use decompInitMod    , only: decompInit_ocn
    use domainMod        , only: domain_check, ldomain, domain_init
    use surfrdMod        , only: surfrd_get_globmask, surfrd_get_grid, surfrd_get_data, surfrd_get_num_patches
    use controlMod       , only: control_init, control_print, NLFilename
    use ncdio_pio        , only: ncd_pio_init
    use initGridCellsMod , only: initGridCells
    use ch4varcon        , only: ch4conrd
    use UrbanParamsType  , only: UrbanInput, IsSimpleBuildTemp
    !
    ! !ARGUMENTS
    integer, intent(in) :: dtime    ! model time step (seconds)

    ! COMPILER_BUG(wjs, 2020-02-20, intel18.0.3) Although gindex_ocn could be
    ! intent(out), intel18.0.3 generates a runtime segmentation fault in runs that don't
    ! have this argument present when this is declared intent(out). (It works fine on
    ! intel 19.0.2 when declared as intent(out).) See also
    ! https://github.com/ESCOMP/CTSM/issues/930.
    integer, pointer, optional, intent(inout) :: gindex_ocn(:)  ! If present, this will hold the decomposition of ocean points (which is needed for the nuopc interface); note that this variable is allocated here, and is assumed to start unallocated
    !
    ! !LOCAL VARIABLES:
    integer           :: ier                     ! error status
    integer           :: i,j,n,k,c,l,g           ! indices
    integer           :: nl                      ! gdc and glo lnd indices
    integer           :: ns, ni, nj              ! global grid sizes
    integer           :: begg, endg              ! processor bounds
    type(bounds_type) :: bounds_proc
    type(bounds_type) :: bounds_clump
    integer           :: nclumps                 ! number of clumps on this processor
    integer           :: nc                      ! clump index
    integer           :: actual_maxsoil_patches  ! value from surface dataset
    integer           :: actual_numcft           ! numcft from sfc dataset
    integer ,pointer  :: amask(:)                ! global land mask
    character(len=32) :: subname = 'initialize1' ! subroutine name
    !-----------------------------------------------------------------------

    call t_startf('clm_init1')

    ! ------------------------------------------------------------------------
    ! Initialize run control variables, timestep
    ! ------------------------------------------------------------------------

    if ( masterproc )then
       write(iulog,*) trim(version)
       write(iulog,*)
       write(iulog,*) 'Attempting to initialize the land model .....'
       write(iulog,*)
       call shr_sys_flush(iulog)
    endif

    call control_init(dtime)
    call ncd_pio_init()
    call surfrd_get_num_patches(fsurdat, actual_maxsoil_patches, actual_numcft)
    call clm_varpar_init(actual_maxsoil_patches, actual_numcft)
    call clm_varcon_init( IsSimpleBuildTemp() )
    call landunit_varcon_init()

    if (masterproc) call control_print()

    call dynSubgridControl_init(NLFilename)

    ! ------------------------------------------------------------------------
    ! Read in global land grid and land mask (amask)- needed to set decomposition
    ! ------------------------------------------------------------------------

    ! global memory for amask is allocate in surfrd_get_glomask - must be
    ! deallocated below
    if (masterproc) then
       write(iulog,*) 'Attempting to read global land mask from ',trim(fatmlndfrc)
       call shr_sys_flush(iulog)
    endif
    call surfrd_get_globmask(filename=fatmlndfrc, mask=amask, ni=ni, nj=nj)

    ! Exit early if no valid land points
    if ( all(amask == 0) )then
       if (masterproc) write(iulog,*) trim(subname)//': no valid land points do NOT run clm'
       noland = .true.
       return
    end if

    ! ------------------------------------------------------------------------
    ! Determine clm gridcell decomposition and processor bounds for gridcells
    ! ------------------------------------------------------------------------

    call decompInit_lnd(ni, nj, amask)
    if (present(gindex_ocn)) then
       call decompInit_ocn(ni, nj, amask, gindex_ocn=gindex_ocn)
    end if
    deallocate(amask)

    if(use_soil_moisture_streams) call decompInit_lnd3D(ni, nj, nlevsoi)
    ! *** Get JUST gridcell processor bounds ***
    ! Remaining bounds (landunits, columns, patches) will be determined
    ! after the call to decompInit_glcp - so get_proc_bounds is called
    ! twice and the gridcell information is just filled in twice

    call get_proc_bounds(begg, endg)

    ! ------------------------------------------------------------------------
    ! Get grid and land fraction (set ldomain)
    ! ------------------------------------------------------------------------

    if (masterproc) then
       write(iulog,*) 'Attempting to read ldomain from ',trim(fatmlndfrc)
       call shr_sys_flush(iulog)
    endif
    call surfrd_get_grid(begg, endg, ldomain, fatmlndfrc)
    if (masterproc) then
       call domain_check(ldomain)
    endif
    ldomain%mask = 1  !!! TODO - is this needed?

    ! Initialize glc behavior
    call glc_behavior%Init(begg, endg, NLFilename)

    ! Initialize urban model input (initialize urbinp data structure)
    ! This needs to be called BEFORE the call to surfrd_get_data since
    ! that will call surfrd_get_special which in turn calls check_urban

    call UrbanInput(begg, endg, mode='initialize')

    ! Allocate surface grid dynamic memory (just gridcell bounds dependent)

    allocate (wt_lunit     (begg:endg, max_lunit           ))
    allocate (urban_valid  (begg:endg                      ))
    allocate (wt_nat_patch (begg:endg, natpft_lb:natpft_ub ))
    allocate (wt_cft       (begg:endg, cft_lb:cft_ub       ))
    allocate (fert_cft     (begg:endg, cft_lb:cft_ub       ))
    allocate (irrig_method (begg:endg, cft_lb:cft_ub       ))
    allocate (wt_glc_mec  (begg:endg, maxpatch_glcmec))
    allocate (topo_glc_mec(begg:endg, maxpatch_glcmec))
    allocate (haslake      (begg:endg                      ))
    ! Read list of Patches and their corresponding parameter values
    ! Independent of model resolution, Needs to stay before surfrd_get_data

    call pftcon%Init()

    ! Read surface dataset and set up subgrid weight arrays
    call surfrd_get_data(begg, endg, ldomain, fsurdat, actual_numcft)

    ! ------------------------------------------------------------------------
    ! Ask Fates to evaluate its own dimensioning needs.
    ! This determines the total amount of space it requires in its largest
    ! dimension.  We are currently calling that the "cohort" dimension, but
    ! it is really a utility dimension that captures the models largest
    ! size need.
    ! Sets:
    ! fates_maxElementsPerPatch
    ! fates_maxElementsPerSite (where a site is roughly equivalent to a column)
    !
    ! (Note: fates_maxELementsPerSite is the critical variable used by CLM
    ! to allocate space)
    ! ------------------------------------------------------------------------

    call set_fates_global_elements(use_fates)

    ! ------------------------------------------------------------------------
    ! Determine decomposition of subgrid scale landunits, columns, patches
    ! ------------------------------------------------------------------------

    call decompInit_clumps(ns, ni, nj, glc_behavior)

    ! *** Get ALL processor bounds - for gridcells, landunit, columns and patches ***

    call get_proc_bounds(bounds_proc)

    ! Allocate memory for subgrid data structures
    ! This is needed here BEFORE the following call to initGridcells
    ! Note that the assumption is made that none of the subgrid initialization
    ! can depend on other elements of the subgrid in the calls below

    call grc%Init  (bounds_proc%begg, bounds_proc%endg)
    call lun%Init  (bounds_proc%begl, bounds_proc%endl)
    call col%Init  (bounds_proc%begc, bounds_proc%endc)
    call patch%Init(bounds_proc%begp, bounds_proc%endp)

    ! Build hierarchy and topological info for derived types
    ! This is needed here for the following call to decompInit_glcp

    call initGridCells(glc_behavior)

    ! Set global seg maps for gridcells, landlunits, columns and patches

    call decompInit_glcp(ns, ni, nj, glc_behavior)

    ! Set filters

    call allocFilters()

    nclumps = get_proc_clumps()
    !$OMP PARALLEL DO PRIVATE (nc, bounds_clump)
    do nc = 1, nclumps
       call get_clump_bounds(nc, bounds_clump)
       call reweight_wrapup(bounds_clump, glc_behavior)
    end do
    !$OMP END PARALLEL DO

    ! ------------------------------------------------------------------------
    ! Remainder of initialization1
    ! ------------------------------------------------------------------------

    ! Set CH4 Model Parameters from namelist.
    ! Need to do before initTimeConst so that it knows whether to
    ! look for several optional parameters on surfdata file.

    if (use_lch4) then
       call ch4conrd()
    end if

    ! Run any requested self-tests
    call self_test_driver(bounds_proc)

    ! Deallocate surface grid dynamic memory for variables that aren't needed elsewhere.
    ! Some things are kept until the end of initialize2; urban_valid is kept through the
    ! end of the run for error checking.

    deallocate (wt_lunit, wt_cft, wt_glc_mec, haslake)

    call t_stopf('clm_init1')

  end subroutine initialize1

  !-----------------------------------------------------------------------
  subroutine initialize2( )
    !
    ! !DESCRIPTION:
    ! CLM initialization - second phase
    !
    ! !USES:

    use shr_orb_mod           , only : shr_orb_decl
    use shr_scam_mod          , only : shr_scam_getCloseLatLon
    use seq_drydep_mod        , only : n_drydep, drydep_method, DD_XLND
    use accumulMod            , only : print_accum_fields
    use clm_varpar            , only : nlevsno
    use clm_varcon            , only : spval
    use clm_varctl            , only : finidat, finidat_interp_source, finidat_interp_dest, fsurdat
    use clm_varctl            , only : use_century_decomp, single_column, scmlat, scmlon, use_cn, use_fates
    use clm_varctl            , only : use_crop, ndep_from_cpl
    use clm_varorb            , only : eccen, mvelpp, lambm0, obliqr
    use clm_time_manager      , only : get_step_size_real, get_curr_calday
    use clm_time_manager      , only : get_curr_date, get_nstep, advance_timestep
    use clm_time_manager      , only : timemgr_init, timemgr_restart_io, timemgr_restart, is_restart
    use CIsoAtmTimeseriesMod  , only : C14_init_BombSpike, use_c14_bombspike, C13_init_TimeSeries, use_c13_timeseries
    use DaylengthMod          , only : InitDaylength
    use dynSubgridDriverMod   , only : dynSubgrid_init
    use dynConsBiogeophysMod  , only : dyn_hwcontent_set_baselines
    use fileutils             , only : getfil
    use initInterpMod         , only : initInterp
    use subgridWeightsMod     , only : init_subgrid_weights_mod
    use histFileMod           , only : hist_htapes_build, htapes_fieldlist, hist_printflds
    use histFileMod           , only : hist_addfld1d, hist_addfld2d, no_snow_normal
    use restFileMod           , only : restFile_getfile, restFile_open, restFile_close
    use restFileMod           , only : restFile_read, restFile_write
    use ndepStreamMod         , only : ndep_init, ndep_interp
    use LakeCon               , only : LakeConInit
    use SatellitePhenologyMod , only : SatellitePhenologyInit, readAnnualVegetation, interpMonthlyVeg
    use SnowSnicarMod         , only : SnowAge_init, SnowOptics_init
    use lnd2atmMod            , only : lnd2atm_minimal
    use NutrientCompetitionFactoryMod, only : create_nutrient_competition_method
    use controlMod            , only : NLFilename
    use clm_instMod           , only : clm_fates
    use BalanceCheckMod       , only : BalanceCheckInit
    !
    ! !ARGUMENTS
    !
    ! !LOCAL VARIABLES:
    integer               :: c,i,j,k,l,p! indices
    integer               :: yr           ! current year (0, ...)
    integer               :: mon          ! current month (1 -> 12)
    integer               :: day          ! current day (1 -> 31)
    integer               :: ncsec        ! current time of day [seconds]
    integer               :: nc           ! clump index
    integer               :: nclumps      ! number of clumps on this processor
    character(len=256)    :: fnamer       ! name of netcdf restart file
    character(len=256)    :: pnamer       ! full pathname of netcdf restart file
    character(len=256)    :: locfn        ! local file name
    type(file_desc_t)     :: ncid         ! netcdf id
    real(r8)              :: dtime        ! time step increment (sec)
    integer               :: nstep        ! model time step
    real(r8)              :: calday       ! calendar day for nstep
    real(r8)              :: caldaym1     ! calendar day for nstep-1
    real(r8)              :: declin       ! solar declination angle in radians for nstep
    real(r8)              :: declinm1     ! solar declination angle in radians for nstep-1
    real(r8)              :: eccf         ! earth orbit eccentricity factor
    type(bounds_type)     :: bounds_proc  ! processor bounds
    type(bounds_type)     :: bounds_clump ! clump bounds
    logical               :: lexist
    integer               :: closelatidx,closelonidx
    real(r8)              :: closelat,closelon
    logical               :: reset_dynbal_baselines_all_columns
    logical               :: reset_dynbal_baselines_lake_columns
    integer               :: begp, endp
    integer               :: begc, endc
    integer               :: begl, endl
    real(r8), pointer     :: data2dptr(:,:) ! temp. pointers for slicing larger arrays
    character(len=32)     :: subname = 'initialize2'
    !----------------------------------------------------------------------

    call t_startf('clm_init2')

    ! ------------------------------------------------------------------------
    ! Determine processor bounds and clumps for this processor
    ! ------------------------------------------------------------------------

    call get_proc_bounds(bounds_proc)
    nclumps = get_proc_clumps()

    ! ------------------------------------------------------------------------
    ! Read in parameters files
    ! ------------------------------------------------------------------------

    call clm_instReadNML( NLFilename )
    allocate(nutrient_competition_method, &
         source=create_nutrient_competition_method(bounds_proc))

    call readParameters(nutrient_competition_method, photosyns_inst)

    ! ------------------------------------------------------------------------
    ! Initialize time manager
    ! ------------------------------------------------------------------------

    if (nsrest == nsrStartup) then
       call timemgr_init()
    else
       call restFile_getfile(file=fnamer, path=pnamer)
       call restFile_open( flag='read', file=fnamer, ncid=ncid )
       call timemgr_restart_io( ncid=ncid, flag='read' )
       call restFile_close( ncid=ncid )
       call timemgr_restart()
    end if

    ! ------------------------------------------------------------------------
    ! Initialize daylength from the previous time step (needed so prev_dayl can be set correctly)
    ! ------------------------------------------------------------------------

    call t_startf('init_orbd')

    calday = get_curr_calday()
    call shr_orb_decl( calday, eccen, mvelpp, lambm0, obliqr, declin, eccf )

    dtime = get_step_size_real()
    caldaym1 = get_curr_calday(offset=-int(dtime))
    call shr_orb_decl( caldaym1, eccen, mvelpp, lambm0, obliqr, declinm1, eccf )

    call t_stopf('init_orbd')

    call InitDaylength(bounds_proc, declin=declin, declinm1=declinm1, obliquity=obliqr)

    ! Initialize Balance checking (after time-manager)
    call BalanceCheckInit()

    ! History file variables

    if (use_cn) then
       call hist_addfld1d (fname='DAYL',  units='s', &
            avgflag='A', long_name='daylength', &
            ptr_gcell=grc%dayl, default='inactive')

       call hist_addfld1d (fname='PREV_DAYL', units='s', &
            avgflag='A', long_name='daylength from previous timestep', &
            ptr_gcell=grc%prev_dayl, default='inactive')
    end if

    ! ------------------------------------------------------------------------
    ! Initialize component data structures
    ! ------------------------------------------------------------------------

    ! Note: new logic is in place that sets all the history fields to spval so
    ! there is no guesswork in the initialization to nans of the allocated variables

    ! First put in history calls for subgrid data structures - these cannot appear in the
    ! module for the subgrid data definition due to circular dependencies that are introduced

    data2dptr => col%dz(:,-nlevsno+1:0)
    col%dz(bounds_proc%begc:bounds_proc%endc,:) = spval
    call hist_addfld2d (fname='SNO_Z', units='m', type2d='levsno',  &
         avgflag='A', long_name='Snow layer thicknesses', &
         ptr_col=data2dptr, no_snow_behavior=no_snow_normal, default='inactive')

    call hist_addfld2d (fname='SNO_Z_ICE', units='m', type2d='levsno',  &
         avgflag='A', long_name='Snow layer thicknesses (ice landunits only)', &
         ptr_col=data2dptr, no_snow_behavior=no_snow_normal, &
         l2g_scale_type='ice', default='inactive')

    col%zii(bounds_proc%begc:bounds_proc%endc) = spval
    call hist_addfld1d (fname='ZII', units='m', &
         avgflag='A', long_name='convective boundary height', &
         ptr_col=col%zii, default='inactive')

    ! If single-column determine closest latitude and longitude

    if (single_column) then
       call getfil (fsurdat, locfn, 0)
       call shr_scam_getCloseLatLon(locfn, scmlat, scmlon, &
            closelat, closelon, closelatidx, closelonidx)
    end if

    ! Initialize instances of all derived types as well as time constant variables
    call clm_instInit(bounds_proc)
    ! Initialize SNICAR optical and aging parameters

    call SnowOptics_init( ) ! SNICAR optical parameters:
    call SnowAge_init( )    ! SNICAR aging   parameters:

    call hist_printflds()

    ! ------------------------------------------------------------------------
    ! Initializate dynamic subgrid weights (for prescribed transient Patches, CNDV
    ! and/or dynamic landunits); note that these will be overwritten in a
    ! restart run
    ! ------------------------------------------------------------------------

    call t_startf('init_dyn_subgrid')
    call init_subgrid_weights_mod(bounds_proc)
    call dynSubgrid_init(bounds_proc, glc_behavior, crop_inst)
    call t_stopf('init_dyn_subgrid')

    ! ------------------------------------------------------------------------
    ! Initialize baseline water and energy states needed for dynamic subgrid operation
    !
    ! This will be overwritten by the restart file, but needs to be done for a cold start
    ! case.
    !
    ! BACKWARDS_COMPATIBILITY(wjs, 2019-03-05) dyn_hwcontent_set_baselines is called again
    ! later in initialization if reset_dynbal_baselines is set. I think we could just have
    ! a single call in that location (adding some logic to also do the call if we're doing
    ! a cold start) once we can assume that all finidat files have the necessary restart
    ! fields on them. But for now, having the extra call here handles the case where the
    ! relevant restart fields are missing from an old finidat file.
    ! ------------------------------------------------------------------------

    !$OMP PARALLEL DO PRIVATE (nc, bounds_clump)
    do nc = 1,nclumps
       call get_clump_bounds(nc, bounds_clump)

       call dyn_hwcontent_set_baselines(bounds_clump, &
            filter_inactive_and_active(nc)%num_icemecc, &
            filter_inactive_and_active(nc)%icemecc, &
            filter_inactive_and_active(nc)%num_lakec, &
            filter_inactive_and_active(nc)%lakec, &
            urbanparams_inst, soilstate_inst, lakestate_inst, water_inst, temperature_inst, &
            reset_all_baselines = .true., &
            ! reset_lake_baselines is irrelevant since reset_all_baselines is true
            reset_lake_baselines = .false.)
    end do
    !$OMP END PARALLEL DO

    ! ------------------------------------------------------------------------
    ! Initialize modules (after time-manager initialization in most cases)
    ! ------------------------------------------------------------------------

    if (use_cn) then
       call bgc_vegetation_inst%Init2(bounds_proc, NLFilename)

       ! NOTE(wjs, 2016-02-23) Maybe the rest of the body of this conditional should also
       ! be moved into bgc_vegetation_inst%Init2

       if (n_drydep > 0 .and. drydep_method == DD_XLND) then
          ! Must do this also when drydeposition is used so that estimates of monthly
          ! differences in LAI can be computed
          call SatellitePhenologyInit(bounds_proc)
       end if

       if ( use_c14 .and. use_c14_bombspike ) then
          call C14_init_BombSpike()
       end if

       if ( use_c13 .and. use_c13_timeseries ) then
          call C13_init_TimeSeries()
       end if
    else
       call SatellitePhenologyInit(bounds_proc)
    end if

    if(use_soil_moisture_streams) then 
       call PrescribedSoilMoistureInit(bounds_proc)
    endif



    ! ------------------------------------------------------------------------
    ! On restart only - process the history namelist.
    ! ------------------------------------------------------------------------

    ! Later the namelist from the restart file will be used.  This allows basic
    ! checking to make sure you didn't try to change the history namelist on restart.

    if (nsrest == nsrContinue ) then
       call htapes_fieldlist()
    end if

    ! ------------------------------------------------------------------------
    ! Read restart/initial info
    ! ------------------------------------------------------------------------

    is_cold_start = .false.
    is_interpolated_start = .false.
    reset_dynbal_baselines_lake_columns = .false.

    if (nsrest == nsrStartup) then

       if (finidat == ' ') then
          if (finidat_interp_source == ' ') then
             is_cold_start = .true.
             if (masterproc) then
                write(iulog,*)'Using cold start initial conditions '
             end if
          else
             if (masterproc) then
                write(iulog,*)'Interpolating initial conditions from ',trim(finidat_interp_source),&
                     ' and creating new initial conditions ', trim(finidat_interp_dest)
             end if
          end if
       else
          if (masterproc) then
             write(iulog,*)'Reading initial conditions from ',trim(finidat)
          end if
          call getfil( finidat, fnamer, 0 )
          call restFile_read(bounds_proc, fnamer, glc_behavior, &
               reset_dynbal_baselines_lake_columns = reset_dynbal_baselines_lake_columns)
       end if

    else if ((nsrest == nsrContinue) .or. (nsrest == nsrBranch)) then

       if (masterproc) then
          write(iulog,*)'Reading restart file ',trim(fnamer)
       end if
       call restFile_read(bounds_proc, fnamer, glc_behavior, &
            reset_dynbal_baselines_lake_columns = reset_dynbal_baselines_lake_columns)
    end if

    ! ------------------------------------------------------------------------
    ! If appropriate, create interpolated initial conditions
    ! ------------------------------------------------------------------------

    if (nsrest == nsrStartup .and. finidat_interp_source /= ' ') then

       is_interpolated_start = .true.

       ! Check that finidat is not cold start - abort if it is
       if (finidat /= ' ') then
          call endrun(msg='ERROR clm_initializeMod: '//&
               'finidat and finidat_interp_source cannot both be non-blank')
       end if

       ! Create new template file using cold start
       call restFile_write(bounds_proc, finidat_interp_dest, writing_finidat_interp_dest_file=.true.)

       ! Interpolate finidat onto new template file
       call getfil( finidat_interp_source, fnamer,  0 )
       call initInterp(filei=fnamer, fileo=finidat_interp_dest, bounds=bounds_proc, &
            glc_behavior=glc_behavior)

       ! Read new interpolated conditions file back in
       call restFile_read(bounds_proc, finidat_interp_dest, glc_behavior, &
            reset_dynbal_baselines_lake_columns = reset_dynbal_baselines_lake_columns)

       ! Reset finidat to now be finidat_interp_dest
       ! (to be compatible with routines still using finidat)
       finidat = trim(finidat_interp_dest)

    end if

    ! ------------------------------------------------------------------------
    ! If requested, reset dynbal baselines
    !
    ! This needs to happen after reading the restart file (including after reading the
    ! interpolated restart file, if applicable).
    ! ------------------------------------------------------------------------

    reset_dynbal_baselines_all_columns = get_reset_dynbal_baselines()
    if (nsrest == nsrBranch) then
       if (reset_dynbal_baselines_all_columns) then
          call endrun(msg='ERROR clm_initializeMod: '//&
               'Cannot set reset_dynbal_baselines in a branch run')
       end if
    else if (nsrest == nsrContinue) then
       ! It's okay for the reset_dynbal_baselines flag to remain set in a continue
       ! run, but we'll ignore it. (This way, the user doesn't have to change their
       ! namelist file for the continue run.)
       reset_dynbal_baselines_all_columns = .false.
    end if
    ! Note that we will still honor reset_dynbal_baselines_lake_columns even in a branch
    ! or continue run: even in these runs, we want to reset those baselines if they are
    ! wrong on the restart file.

    if (masterproc) then
       if (reset_dynbal_baselines_all_columns) then
          write(iulog,*) ' '
          write(iulog,*) 'Resetting dynbal baselines for all columns'
          write(iulog,*) ' '
       else if (reset_dynbal_baselines_lake_columns) then
          write(iulog,*) ' '
          write(iulog,*) 'Resetting dynbal baselines for lake columns'
          write(iulog,*) ' '
       end if
    end if

    !$OMP PARALLEL DO PRIVATE (nc, bounds_clump)
    do nc = 1,nclumps
       call get_clump_bounds(nc, bounds_clump)

       call dyn_hwcontent_set_baselines(bounds_clump, &
            filter_inactive_and_active(nc)%num_icemecc, &
            filter_inactive_and_active(nc)%icemecc, &
            filter_inactive_and_active(nc)%num_lakec, &
            filter_inactive_and_active(nc)%lakec, &
            urbanparams_inst, soilstate_inst, lakestate_inst, &
            water_inst, temperature_inst, &
            reset_all_baselines = reset_dynbal_baselines_all_columns, &
            reset_lake_baselines = reset_dynbal_baselines_lake_columns)
    end do
    !$OMP END PARALLEL DO

    ! ------------------------------------------------------------------------
    ! Initialize nitrogen deposition
    ! ------------------------------------------------------------------------

    if (use_cn) then
       call t_startf('init_ndep')
       if (.not. ndep_from_cpl) then
          call ndep_init(bounds_proc, NLFilename)
          call ndep_interp(bounds_proc, atm2lnd_inst)
       end if
       call t_stopf('init_ndep')
    end if

    ! ------------------------------------------------------------------------
    ! Initialize active history fields.
    ! ------------------------------------------------------------------------

    ! This is only done if not a restart run. If a restart run, then this
    ! information has already been obtained from the restart data read above.
    ! Note that routine hist_htapes_build needs time manager information,
    ! so this call must be made after the restart information has been read.

    if (nsrest /= nsrContinue) then
       call hist_htapes_build()
    end if

    ! ------------------------------------------------------------------------
    ! Initialize variables that are associated with accumulated fields.
    ! ------------------------------------------------------------------------

    ! The following is called for both initial and restart runs and must
    ! must be called after the restart file is read

    call atm2lnd_inst%initAccVars(bounds_proc)
    call temperature_inst%initAccVars(bounds_proc)
    call water_inst%initAccVars(bounds_proc)
    call energyflux_inst%initAccVars(bounds_proc)
    call canopystate_inst%initAccVars(bounds_proc)

    call bgc_vegetation_inst%initAccVars(bounds_proc)

    if (use_crop) then
       call crop_inst%initAccVars(bounds_proc)
    end if

    !------------------------------------------------------------
    ! Read monthly vegetation
    !------------------------------------------------------------

    ! Even if CN is on, and dry-deposition is active, read CLMSP annual vegetation
    ! to get estimates of monthly LAI

    if ( n_drydep > 0 .and. drydep_method == DD_XLND )then
       call readAnnualVegetation(bounds_proc, canopystate_inst)
       if (nsrest == nsrStartup .and. finidat /= ' ') then
          ! Call interpMonthlyVeg for dry-deposition so that mlaidiff will be calculated
          ! This needs to be done even if CN or CNDV is on!
          call interpMonthlyVeg(bounds_proc, canopystate_inst)
       end if
    end if

    !------------------------------------------------------------
    ! Determine gridcell averaged properties to send to atm
    !------------------------------------------------------------

    if (nsrest == nsrStartup) then
       call t_startf('init_map2gc')
       call lnd2atm_minimal(bounds_proc, &
            water_inst, surfalb_inst, energyflux_inst, lnd2atm_inst)
       call t_stopf('init_map2gc')
    end if

    !------------------------------------------------------------
    ! Initialize sno export state to send to glc
    !------------------------------------------------------------

    !$OMP PARALLEL DO PRIVATE (nc, bounds_clump)
    do nc = 1,nclumps
       call get_clump_bounds(nc, bounds_clump)

       call t_startf('init_lnd2glc')
       call lnd2glc_inst%update_lnd2glc(bounds_clump,       &
            filter(nc)%num_do_smb_c, filter(nc)%do_smb_c,   &
            temperature_inst, water_inst%waterfluxbulk_inst, topo_inst, &
            init=.true.)
       call t_stopf('init_lnd2glc')
    end do
    !$OMP END PARALLEL DO

    !------------------------------------------------------------
    ! Deallocate wt_nat_patch
    !------------------------------------------------------------

    ! wt_nat_patch was allocated in initialize1, but needed to be kept around through
    ! initialize2 for some consistency checking; now it can be deallocated

    deallocate(wt_nat_patch)

    ! --------------------------------------------------------------
    ! Initialise the fates model state structure
    ! --------------------------------------------------------------

    if ( use_fates .and. .not.is_restart() .and. finidat == ' ') then
       call clm_fates%init_coldstart(water_inst%waterstatebulk_inst, &
            water_inst%waterdiagnosticbulk_inst, canopystate_inst, &
            soilstate_inst)
    end if

    ! topo_glc_mec was allocated in initialize1, but needed to be kept around through
    ! initialize2 because it is used to initialize other variables; now it can be
    ! deallocated

    deallocate(topo_glc_mec, fert_cft, irrig_method)

    !------------------------------------------------------------
    ! Write log output for end of initialization
    !------------------------------------------------------------

    call t_startf('init_wlog')
    if (masterproc) then
       write(iulog,*) 'Successfully initialized the land model'
       if (nsrest == nsrStartup) then
          write(iulog,*) 'begin initial run at: '
       else
          write(iulog,*) 'begin continuation run at:'
       end if
       call get_curr_date(yr, mon, day, ncsec)
       write(iulog,*) '   nstep= ',get_nstep(), ' year= ',yr,' month= ',mon,&
            ' day= ',day,' seconds= ',ncsec
       write(iulog,*)
       write(iulog,'(72a1)') ("*",i=1,60)
       write(iulog,*)
    endif
    call t_stopf('init_wlog')

    call t_stopf('clm_init2')

    if (water_inst%DoConsistencyCheck()) then
       !$OMP PARALLEL DO PRIVATE (nc, bounds_clump)
       do nc = 1,nclumps
          call get_clump_bounds(nc, bounds_clump)
          call water_inst%TracerConsistencyCheck(bounds_clump, 'end of initialization')
       end do
       !$OMP END PARALLEL DO
    end if

  end subroutine initialize2

end module clm_initializeMod
